//#ist=0x8088cac) at /usr/include/c++/4.6/bits/stl_construct.h:127
#include "CSpaceCharge.h"
#include <iostream>
#include <cmath>
#include <stdio.h>
#include <iomanip>
#include <sstream>
using namespace std;
const double dzero=1e-9;
const double dPermit=8.85e-12;
const double dElecCharge=1.602e-19;
//#include "common.h"

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])
//电流单位是mA
CSpaceCharge::CSpaceCharge()
{}
CSpaceCharge::~CSpaceCharge()
{
  delete dEnergyK;
}
void CSpaceCharge::setFlag(int _iflag)
{
  iflag=_iflag;
}
void CSpaceCharge::formGrid()
{
  /*
     plim(&dPhis, &dpLeft);
     double pright=-dPhis;
     double cay=2*pi/dLCell;
     double zmax=(pright-dpLeft)/cay;  
     double ddz=zmax/nz;
  //    if(dz<ddz)dz=ddz;    //dz is the long grid length,while ddz is the value after regulate the grid
   */
  vvdaa.resize(nr);
  vvder.resize(nr);
  vvdez.resize(nr);
  compData.resize(nr);
  vvders.resize(nr*nr);
  vvdezs.resize(nr*nr);
  for(int i=0;i<vvders.size();i++)
  {
    vvders[i].resize(nz);
    vvdezs[i].resize(nz);
    for(int j=0;j<vvders[i].size();j++)
    {
      vvders[i][j]=0;
      vvdezs[i][j]=0;
    }
  }
  for(int i=0;i<vvdaa.size();i++)
  {
    vvdaa[i].resize(nz);
    vvder[i].resize(nz);
    vvdez[i].resize(nz);
    compData[i].resize(nz*2);
    for(int k=0;k<vvdaa[i].size();k++)
    {
      vvdaa[i][k]=0;
      vvder[i][k]=0;
      vvdez[i][k]=0;
      compData[i][2*k]=0;
      compData[i][2*k+1]=0;
    }
  }  

  vdrm.resize(nr);
  vdzm.resize(nz);
  vdrs.resize(nr);
  vdzs.resize(nz);
  rHalf=0.0;
  zHalf=(nz-1)*dz/2.0;
  if(iflag==1)
  {
    rHalf=(nr-1)*dr/2.0;
  }
  if(iflag==1)
  {
    for(int i=0;i<nr;i++)
    {
      vdrs[i]=i*dr;
      vdrm[i]=vdrs[i]+dr/2.0;
    }
  }
  else if(iflag==0)
  {
    vdrm[0]=0;
    for(int i=1;i<nr;i++)
    {
      vdrm[i]=i*dr;
      vdrs[i-1]=sqrt(0.5*(pow(vdrm[i-1],2)+pow(vdrm[i],2)));
    }
    vdrs[nr-1]=vdrm[nr-1]+dr/2.0;
  }
  for(int i=0;i<nz;i++)
  {
    vdzm[i]=i*dz;
    vdzs[i]=vdzm[i]+dz/2.0;
  }
}
void CSpaceCharge::initField()       //初始化电场
{
  double *er1, *ez1;
  //c1=572167*dQ/dErest/1e6*(dFreq/1e6*dWaveLength*100)*(dFreq/1e6*dWaveLength*100);
  //c1单位是cm3/um2
  //c1=572167*dQ/dErest*1e-6*clight*clight/1e8; //clight的单位换算成是cm/us.和上面的结果应该是一样的  
  c1=dQ/2.0/pi/pi/dPermit;
  //(2/pi)*e/(4*pi*epsilon)=572167 cm mev/coul
  er1 = new double;
  ez1 = new double;
  for (int k=0; k<nr; k++)
  {
    for(int i=0; i<nr; i++)
    {
      for(int j=0; j<nz; j++)
      {
        flds(vdrs[k]/100.0, dz/2.0/100.0, nip, vdrm[i]/100.0, vdzm[j]/100.0, pl, er1, ez1);
        vvders[k*nr+i][j] = c1*(*er1);
        vvdezs[k*nr+i][j] = c1*(*ez1);
      }
    }
  }
  delete er1;
  delete ez1;
}
void CSpaceCharge::scheff(double dist)
{
  dt=dist;
  dQ =dI/(dFreq*iParticleNumTotal);  // I is average current in one rf period
  c1=dQ/2.0/pi/pi/dPermit;
  //c1=572167*dQ/dErest/1e6*(dFreq/1e6*dWaveLength*100)*(dFreq/1e6*dWaveLength*100);
  //c1单位是cm3/um2
  //c1=572167*dQ/dErest*1e-6*clight*clight/1e8; //clight的单位是cm/us.和上面的结果应该是一样的  
  //c2=dWaveLength*100/2/pi;
  //c3=dt;	// schefftm(htau); dt=htau is time length and the utility is us（微秒）.
  //c3=dLCell;
  /*
     if(fabs(dt)<dzero)
     {
  //	    dz = longigridlength;
  //	double pleft;
  formGrid();
  }
   */
  //else
  if(iflag==0) initField();       	
  /*
  ofstream ofr("vvders");
  ofstream ofz("vvdezs");
  for(int i=0;i<nz;i++)
  {
    ofr<<"particle "<<i<<endl;
    ofz<<"particle "<<i<<endl;
    for(int j=0;j<nr;j++)
    {
      for(int k=0;k<nz;k++)
      {
        if(k>i)
        {
          ofr<<vvders[j][abs(k-i)]<<" ";
          ofz<<vvdezs[j][abs(k-i)]<<" ";
          }
        else
        {
          ofr<<vvders[j][abs(k-i)+1]<<" ";
          ofz<<-vvdezs[j][abs(k-i)+1]<<" ";
          }

      }
      ofr<<endl;
      ofz<<endl;
    }
    ofr<<endl;
    ofz<<endl;
  }
  */
  //getvvd(vvders "vvders",0,1);
  //getvvd(vvders "vvderz",0,1);
  double xbar=0;	//MainCode.h
  double ybar=0;	//MainCode.h
  double xsq=0;
  double ysq=0;
  double pbar=0;	//pbar 为正离子纵向坐标z的平均值	//MainCode.h
  double pmax=0;
  double pmin=0;

  for(int np=0; np<cord->size(); np++)
  {
    xbar=xbar+(*cord)[np][0];
    ybar=ybar+(*cord)[np][2];
    xsq=xsq+(*cord)[np][0]*(*cord)[np][0];
    ysq=ysq+(*cord)[np][2]*(*cord)[np][2];
    pbar=pbar+(*cord)[np][4];
    if((*cord)[np][4]<pmin) pmin=(*cord)[np][4];
    if((*cord)[np][4]>pmax) pmax=(*cord)[np][4];
  }  
  /*
     double plT=(pmax-pmin);
     double zlT=(pmax-pmin)/2/pi*dWaveLength*dBetaS;
     if(dz*nz>3*zlT) 
     {dpLeft=pmin-plT;dz=3*zlT/nz; 
     for(int i=0; i<nz+1; i++)
     {
     vdzm[i]=i*dz;
     vdzzs[i]=vdzm[i]+dz/2;
     }
     }
   */
  ng=cord->size();                                     //Rx---Right-----i
  xbar=xbar/ng;                                        //Lx---Left------i-1
  ybar=ybar/ng;                                        //Uy---Up--------j
  xsq=xsq/ng;                                          //Dy---Down------j-1
  ysq=ysq/ng;                                          //Fz---forward---k 
  pbar=pbar/ng;                                        //Bz---BackWard--k-1

  double dxP,dyP,dzP,dphiP,drP;
  int    ir1,ir2,iz1,iz2;
  double dr1,dr2,dz1,dz2;
  int    iSign=1;
//  ofstream ofr("beamTrans");
  for(int np=0;np<cord->size();np++)
  {
    dxP=(*cord)[np][0]; dyP=(*cord)[np][2];dphiP=(*cord)[np][4];
    if(dxP<0&&iflag==1)  iSign=-1;
    dzP=dphiP*dBetaS*dWaveLength/2.0/pi+zHalf;
    //drP=sqrt(dxP*dxP+dyP*dyP)*iSign;
    drP=sqrt(dxP*dxP+dyP*dyP)*iSign+rHalf;
    iSign=1;
    ir1=drP/dr;iz1=dzP/dz;
    ir2=ir1+1;iz2=iz1+1;
  //  ofr<<drP-rHalf<<" "<<dzP-zHalf<<endl;
    if(iflag==0)
    {
      if(drP<vdrs[ir1]&&ir1!=0) ir2=ir1-1;
      else if(drP<vdrs[ir1]&&ir1==0)  {ir1=0;ir2=0;}
      if(dzP<vdzs[iz1]&&iz1!=0) iz2=iz1-1;
      else if(dzP<vdzs[iz1]&&iz1==0) {iz2=0;iz1=0;};
    }
    if(ir1<0) {ir1=0;ir2=0;}
    if(ir1>=nr-1){ir1=nr-1;ir2=nr-1;}
    if(iz1<0) {iz1=0;iz2=0;}
    if(iz1>=nz-1){iz1=nz-1;iz2=nz-1;}

    dr2=1.0;dz2=1.0;
    if(ir1!=ir2)
    {
      dr2=(drP-vdrs[ir1])/(vdrs[ir2]-vdrs[ir1]);
    }
    if(iz1!=iz2)
    {
      dz2=(dzP-vdzs[iz1])/(vdzs[iz2]-vdzs[iz1]);
    }
    dr1=1.0-dr2;
    dz1=1.0-dz2;
    /*
    vvdaa[ir1][iz1]=vvdaa[ir1][iz1]+dr2*dz2;
    vvdaa[ir1][iz2]=vvdaa[ir1][iz2]+dr2*dz1;
    vvdaa[ir2][iz1]=vvdaa[ir2][iz1]+dr1*dz2;
    vvdaa[ir2][iz2]=vvdaa[ir2][iz2]+dr1*dz1;
    */
    vvdaa[ir1][iz1]=vvdaa[ir1][iz1]+dr1*dz1;
    vvdaa[ir1][iz2]=vvdaa[ir1][iz2]+dr1*dz2;
    vvdaa[ir2][iz1]=vvdaa[ir2][iz1]+dr2*dz1;
    vvdaa[ir2][iz2]=vvdaa[ir2][iz2]+dr2*dz2;
  }
  //test===============================
 // ofr.close();
  getvvd(vvdaa,"vvdaa",0,0);
  /*
  ofstream of("vvdaa");
  for(int i=0;i<vvdaa.size();i++)
  {
    for(int j=0;j<vvdaa[i].size();j++)
      of<<vvdaa[i][j]<<" ";
    of<<endl;
  }
  of.close();
  */
  //end test===========================

  for(int i=0;i<vvdaa.size();i++)
    for(int j=0;j<vvdaa[i].size();j++)
    {
      compData[i][2*j]=vvdaa[i][j]*dQ/dr/dz*1.0e4;    //FFT2D???????????????????????
      compData[i][2*j+1]=0;    //FFT2D???????????????????????
    }
//  getPotential("comp-after-aa");
  if(iflag==0) EFieldMap();
  if(iflag==1) EFieldMap2();
}

void CSpaceCharge::Boundary(int _iflag)
{
  vismax.resize(nz);			
  viemax.resize(nz);			
  for(int i=0;i<nz;i++)
  {
    vismax[i]=nr-1;
    viemax[i]=nr-1;
  }
  /*
     int i,j;
     if(_iflag==0)
     {
     for(j=0; j<nz; j++)	
     {
     for(i=nr-1; i>0;)
     {
     if(vvdaa[i][j]==0) i=i-1;
     if(vvdaa[i][j]!=0) break;
     }
     vismax[j]=i+1;  //i是角标的编号，其实相当于是第i+1个格点
     }
     viemax[0]=vismax[0];
     for(j=1; j<nz-1; j++)
     viemax[j]=1+max(vismax[j+1],vismax[j]);
     viemax[nz-1]=vismax[nz-1];
     }
   */
}
void CSpaceCharge::EFieldMap()      //求解每个场点网格点电场值
{
  Boundary(0);
  int js1,ism,iem,iTmp;
  double aa;// aa表示当前格点电量加权值
  for(int js=0; js<nz; js++)		//正离子
  {
    ism=vismax[js];
    //		ism 为纵向第 js 层网格中的电荷环数
    if(ism==0) continue;
    for(int is=0; is<ism; is++)
    {
      aa=vvdaa[is][js];
//        aa=1.0;
      if (aa==0)continue;
      for(int je=0;je<nz;je++)
      {
        iem=viemax[je];
        for(int ie=0;ie<iem;ie++)
        {
          if(je>js)
          {
            vvdez[ie][je]=vvdez[ie][je]+aa*vvdezs[is*nr+ie][abs(je-js)];
            vvder[ie][je]=vvder[ie][je]+aa*vvders[is*nr+ie][abs(je-js)];
          }
          else
          {
            vvder[ie][je]=vvder[ie][je]+aa*vvders[is*nr+ie][abs(je-js)+1];
            vvdez[ie][je]=vvdez[ie][je]-aa*vvdezs[is*nr+ie][abs(je-js)+1];
          }
        }
      }
    }
  }
  getvvd(vvder,"vvder",0,1);
  getvvd(vvdez,"vvdez",0,1);
  ofstream ofe("vvde");
  for(int i=0;i<vvder.size();i++)
    for(int j=0;j<vvder[i].size();j++)
    {
      ofe<<vdrm[i]<<" "<<vdzm[j]<<" "<<sqrt(vvder[i][j]*vvder[i][j]+vvdez[i][j]*vvdez[i][j])<<endl;
    }
  ofe.close();
}

void CSpaceCharge::EFieldMap2()      //求解每个场点网格点电场值
{
  //FFT2D(&compData);
  FFT1D(&compData);
  //getPotential("pot-after-FFT2D");
  double dAngle;
  for(int i=0;i<compData.size();i++)
  {
    for(int k=0;k<compData[i].size();k=k+2)
    {
      dAngle=atan2(compData[i][k+1],compData[i][k]);
      compData[i][k]=sqrt(compData[i][k]*compData[i][k]+compData[i][k+1]*compData[i][k+1]);
      compData[i][k+1]=dAngle;
    }
  }
  Potential();              
  //getPotential("pot-after-pot");
  double dAmp;
  for(int i=0;i<compData.size();i++)
  {
    for(int j=0;j<compData[i].size();j=j+2)
    {
      dAngle=compData[i][j+1];
      dAmp=compData[i][j];
      compData[i][j+1]=dAmp*sin(dAngle);
      compData[i][j]=dAmp*cos(dAngle);
    }
  }
  IFFT1D(&compData);
  //IFFT2D(&compData);
  /*
  for(int i=0;i<compData.size();i++)
  {
    for(int k=0;k<compData[i].size();k=k+2)
    {
      dAngle=atan2(compData[i][k+1],compData[i][k]);
      compData[i][k]=sqrt(compData[i][k]*compData[i][k]+compData[i][k+1]*compData[i][k+1]);
      compData[i][k+1]=dAngle;
    }
  }
  */
  getPotential("pot");
  double dBU,dBD,dFU,dFD;
  for(int i=0;i<vdrm.size()-1;i++)
    for(int k=0;k<vdzm.size()-1;k++)
    {
      dBU=compData[i+1][2*k];
      dBD=compData[i][2*k];
      dFU=compData[i+1][2*(1+k)];
      dFD=compData[i][2*(k+1)];
      vvder[i][k]=(dBU+dFU-dBD-dFD)/2.0/dr*100;
      vvdez[i][k]=(dFU+dFD-dBU-dBD)/2.0/dz*100;
    }
    getvvd(vvder,"vvder-FFT",0,1);
    getvvd(vvdez,"vvdez-FFT",0,1);
}
//getvvdez();
//getvvder();
//after IFFT2D,the useful part is the real part. This can be understood from amplitude and angle mode.
//where the angle mean a "zhi hou".
//复数的电荷密度仅仅实部有用!!!!!!!!!!
//下面是把单位和flag=0时候统一起来

void CSpaceCharge::FFT2D(vector<vector<double> > *vvdComp)
{
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double aTmp[2*nz];
  for(int j=0;j<nr;j++)
  {
    for(int k=0;k<nz;k++)
    {
      REAL(aTmp,k)=(*vvdComp)[j][2*k];
      IMAG(aTmp,k)=(*vvdComp)[j][2*k+1];
    }
    gsl_fft_complex_forward(aTmp,1,nz,compwave,compwork);
    for(int k=0;k<nz;k++)
    {
      (*vvdComp)[j][2*k]=REAL(aTmp,k);
      (*vvdComp)[j][2*k+1]=IMAG(aTmp,k);
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
  compwork=gsl_fft_complex_workspace_alloc(nr);
  compwave=gsl_fft_complex_wavetable_alloc(nr);
  double bTmp[2*nr];
  for(int i=0;i<nz;i++)
  {
    for(int k=0;k<nr;k++)
    {
      REAL(bTmp,k)=(*vvdComp)[k][2*i];
      IMAG(bTmp,k)=(*vvdComp)[k][2*i+1];
    }
    gsl_fft_complex_forward(bTmp,1,nr,compwave,compwork);
    for(int k=0;k<nr;k++)
    {
      //      (*vvdComp)[j][i]=bTmp[2*j];
      //      (*vvdComp)[j][i+1]=bTmp[2*j+1];
      (*vvdComp)[k][2*i]=REAL(bTmp,k);
      (*vvdComp)[k][2*i+1]=IMAG(bTmp,k);
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
}
void CSpaceCharge::FFT1D(vector<vector<double> > *vvdComp)
{
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double aTmp[2*nz];
  for(int i=0;i<nr;i++)
  {              
    for(int k=0;k<nz;k++)
    {
      REAL(aTmp,k)=(*vvdComp)[i][2*k];
      IMAG(aTmp,k)=(*vvdComp)[i][2*k+1];
    }
    gsl_fft_complex_forward(aTmp,1,nz,compwave,compwork);
    //    for(int k=0;k<2*(nz+1);k=k+2)
    for(int k=0;k<nz;k++)
    {
      (*vvdComp)[i][2*k]=REAL(aTmp,k);
      (*vvdComp)[i][2*k+1]=IMAG(aTmp,k);
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
}
void CSpaceCharge::IFFT2D(vector<vector<double> >* vvdComp)
{
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork=gsl_fft_complex_workspace_alloc(nr);
  compwave=gsl_fft_complex_wavetable_alloc(nr);
  double aTmp[2*nr];
  for(int i=0;i<nz;i++)
  {
    for(int k=0;k<nr;k++)
    {
      REAL(aTmp,k)=(*vvdComp)[k][2*i];
      IMAG(aTmp,k)=(*vvdComp)[k][2*i+1];
    }
    gsl_fft_complex_inverse(aTmp,1,nr,compwave,compwork);
    for(int k=0;k<nr;k++)
    {
      (*vvdComp)[k][2*i]=REAL(aTmp,k);
      (*vvdComp)[k][2*i+1]=IMAG(aTmp,k);
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);

  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double cTmp[2*nz];
  for(int i=0;i<nr;i++)
  {
    for(int k=0;k<nz;k++)
    {
      REAL(cTmp,k)=(*vvdComp)[i][2*k];
      IMAG(cTmp,k)=(*vvdComp)[i][2*k+1];
    }
    gsl_fft_complex_inverse(cTmp,1,nz,compwave,compwork);
    for(int k=0;k<nz;k++)
    {
      (*vvdComp)[i][2*k]=REAL(cTmp,k);
      (*vvdComp)[i][2*k+1]=IMAG(cTmp,k);
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
}
void CSpaceCharge::IFFT1D(vector<vector<double> >* vvdComp)
{
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double cTmp[2*nz];
  for(int i=0;i<nr;i++)
  {
    for(int k=0;k<nz;k=k++)
    {
      REAL(cTmp,k)=(*vvdComp)[i][2*k];
      IMAG(cTmp,k)=(*vvdComp)[i][2*k+1];
    }
    gsl_fft_complex_inverse(cTmp,1,nz,compwave,compwork);
    for(int k=0;k<nz;k++)
    {
      (*vvdComp)[i][2*k]=REAL(cTmp,k);
      (*vvdComp)[i][2*k+1]=IMAG(cTmp,k);
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
}
void CSpaceCharge::Potential()       
{     
  double dkr[nr],dkz[nz],dif2K[nr][nz];
  double dTmpr,dTmpz;
  double dif;
  dkr[0]=0;dkz[0]=0;
  //FFT2D============================================
  /*
  for(int i=1;i<nr;i++)
  {
    if(i<nr/2)
      dkr[i]=i*(2*pi/dr/nr);
    else dkr[i]=(i-nr)*(2*pi/dr/nr);
    dTmpr=sin(dkr[i]*dr)/(dkr[i]*dr);

    for(int j=1;j<nz;j++)
    {
      if(j<nz/2)
        dkz[j]=j*(2*pi/dz/nz);
      else       dkz[j]=(j-nz)*(2*pi/dz/nz);
      dTmpz=sin(dkz[j]*dz)/(dkz[j]*dz);
      dif=dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4+dTmpr*dTmpr*dkr[i]*dkr[i]*1.0e4;
      compData[i][2*j]=compData[i][2*j]/dif/dPermit;
    }
  }
  */
  //FFT1D============================================
  /*
  for(int j=0;j<nz;j++)
  {
    if(j<nz/2)
           //dkz[j]=j*(2*pi/10.0);
      dkz[j]=(1+j)*(2*pi/dz/nz);
    //else       dkz[j]=(j-nz)*(2*pi/10.0);
    else       dkz[j]=(j-nz)*(2*pi/dz/nz);
    dTmpz=sin(dkz[j]*dz)/(dkz[j]*dz);
    dif=dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4;
    for(int i=0;i<nr;i++)
      compData[i][2*j]=compData[i][2*j]/dif/dPermit;
  }
  */
  for(int j=1;j<nz;j++)
  {
    if(j<nz/2+1)
           //dkz[j]=j*(2*pi/10.0);
      dkz[j]=j*(2*pi/dz/nz);
    //else       dkz[j]=(j-nz)*(2*pi/10.0);
    else       dkz[j]=(j-nz)*(2*pi/dz/nz);
    dTmpz=sin(dkz[j]*dz)/(dkz[j]*dz);
    dif=dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4;
    for(int i=1;i<nr;i++)
      compData[i][2*j]=compData[i][2*j]/dif/dPermit;
      compData[0][2*j]=0;
  }
    for(int i=0;i<nr;i++)
      compData[i][0]=0;
  /*
  for(int i=0;i<nr;i++)
  {
    for(int j=0;j<nz;j++)
      compData[i][2*j]=compData[i][2*j]/dif/dPermit;
//    compData[i][0]=0;
  }
  */
}
void CSpaceCharge::setFreq(double _dFreq)
{
  dFreq=_dFreq*1e6;
  dWaveLength=clight/dFreq;
}
void CSpaceCharge::InputBeam(vector<vector<double> >* _vvdBeam,double _dI)
{
  cord=_vvdBeam;                                                    
  dI=_dI/1000.0;
  iParticleNumTotal=_vvdBeam->size();      //第一个粒子是同步粒子，不包括在粒子数范围内
  dEnergyK=new double[iParticleNumTotal];
  for(int n=0;n<iParticleNumTotal;n++)
    dEnergyK[n]=(*_vvdBeam)[n][5];
  double dEnergyS=dEnergyK[0];
  double dTmp=dEnergyS/dErest+1;
  //    double dTmp=dEnergyS/dErest/1e6+1;
  dBetaS=sqrt(1-1/dTmp/dTmp);
  dLCell=dBetaS*dWaveLength;
}
void CSpaceCharge::InGrid(int _nr,int _nz,double _dr,double _dz)
{
  dr=_dr;
  dz=_dz;
  nz=_nz;
  nr=_nr;
  if(iflag==1)
  {
    nr=nr*2;           //in cm
    nz=nz*2;
  }
  /*
     dDensityK=new double*[nr];
     dDensity=new double*[nr];
     for(int iTmp=0;iTmp<nr;iTmp++)
     {
     dDensityK[iTmp]=new double[nz+1];
     dDensity[iTmp]=new double[nz+1];
     }
   */
  /*
     compData=new double*[nr+1];
     for(int i=0;i<nr;i++)
     compData[i]=new double[2*(nz+1)];
   */
  /*
     compData.resize(nr);
     for(int i=0;i<compData.size();i++)
     compData[i].resize(2*(nz+1));
   */
}
void CSpaceCharge::flds(double r, double z, int nip, double r1, double z1, double pl, double* er, double* ez)
{
  double d=z1-z;
  double c=pow((r1-r),2);
  double b=pow((r1+r),2);
  double a=4.0*r*r1/(b+d*d);
  double *ee, *ek;
  ee=new double;
  ek=new double;
  eint(a, ee, ek);
  double er1=0.0;
  double ez1=0.0;
  a=sqrt(b+d*d);

  if (r1!=0)
    er1=((*ek)-(r*r-r1*r1+d*d)*(*ee)/(c+d*d))/(2.0*r1*a);
  ez1=d*(*ee)/(a*(c+d*d));
  if (nip!=0)
  {
    for(int i=0; i<nip; i++)
    {
      int xi=i+1;
      for(int j=0; j<2; j++)
      {
        d=z1-(z+xi*pl);

        if(fabs(d)<dz/2.0)
        {
          if(d>=0)d=dz/2.0;
          if(d<0)d=-1.0*dz/2.0;
        }

        a=4.0*r*r1/(b+d*d);
        eint (a,ee,ek);
        a=sqrt(b+d*d);
        if (r1!=0)
          er1=er1+((*ek)-(r*r-r1*r1+d*d)*(*ee)/(c+d*d))/(2.0*r1*a);
        ez1=ez1+d*(*ee)/(a*(c+d*d));
        xi=-xi;
      }
    }
  }
  *er=er1;
  *ez=ez1;
  delete ee;
  delete ek;
}
void CSpaceCharge::eint(double a, double *ee, double *ek)
{

  double b=1.0-a;
  double c=log(b);
  *ee=1.0+b*(0.4630106-0.2452740*c+
      b*(0.1077857-0.04125321*c));
  *ek=1.38629436-0.5*c+b*(0.1119697-0.1213486*c+
      b*(0.07253230-0.0288747*c));
}
void CSpaceCharge::setPhis(double _dPhis)   //attention in degree
{
  dPhis=_dPhis/180*pi;
}
void CSpaceCharge::InputStruct(double _dFreq,double _dErest)
{
  dFreq=_dFreq*1.0e6;
  dWaveLength=clight/dFreq;
  dErest=_dErest;
}
/*
   void CSpaceCharge::Interpret()
   {
   zpTozt(cord,dFreq);
   double dGam;
   double drP,dzP;
   double dLeft,dRight,dLow,dHigh;
   double drb,dzb;
   int  iHigh,iRight,iLeft,iLow;
   int l,m;
   ofstream of("czd");
   of<<" iHigh "<<" iRight "<<" dHigh "<<"dLow"<<" dLeft"<<" dRight"<<" vvdez[iLow][iLeft] "<<"vvdez[iLow][iRight]"<<" vvdez[iHigh][iLeft] "<<" vvdez[iHigh][iRight] "<<" czd "<<endl;
   for(np=1; np!=cord->size(); np++)
   {
   dGam=dEnergyK[np]/dErest+1.0;
   drP=sqrt(pow((*cord)[np][0]-xbar,2)+pow((*cord)[np][2]-ybar,2));	
   if(drP==0)drP=.000001;
   xOverR=((*cord)[np][0]-xbar)/drP;
   yOverR=((*cord)[np][2]-ybar)/drP;
   dzP=((*cord)[np][4]-dpLeft);
   dzP=dzP/2/pi*dBetaS*dWaveLength;
//if((drP<=rmax)&&(fabs(dzP)<=dz*nz))	//求所有离子感受到的正离子的力
if((drP<=rmax))	//求所有离子感受到的正离子的力
{
drb=drP/dr;
iLow=drb;
iHigh=1+drb;
dLow=drb-iLow;    //a
dHigh=1.0-dLow;       //b

//      dzb=(dzP+nz*dz/2)/dz;
dzb=dzP/dz;
iRight=1+dzb;
iLeft=iRight-1;
dLeft=dzb-(iRight-1);  //c
dRight=1.0-dLeft;  //d

//      l=iHigh+(iRight-1)*(1+nr);           //depends!!!!!
//      m=l+(nr+1);
//c3=dt means 计算的空间电荷作用时间
//crp means c3*格点电场/gamma3，means rp的变化量
//czd means c3*格点电场/gamma3. 表示动能的变化量
//c3的单位是us
crp=c3*(dRight*(dLow*vvder[iHigh][iLeft]+dHigh*vvder[iLow][iLeft])+dLeft*(dLow*vvder[iHigh][iRight]+dHigh*vvder[iLow][iRight]))/dGam/dGam/dGam;
czd=c3*(dRight*(dLow*vvdez[iHigh][iLeft]+dHigh*vvdez[iLow][iLeft])+dLeft*(dLow*vvdez[iHigh][iRight]+dHigh*vvdez[iLow][iRight]))/dGam;
of<<iHigh<<" "<<iRight<<" "<<dHigh<<" "<<dLow<<" "<<dLeft<<" "<<dRight<<" "<<vvdez[iLow][iLeft]<<" "<<vvdez[iLow][iRight]<<" "<<vvdez[iHigh][iLeft]<<" "<<vvdez[iHigh][iRight]<<" "<<czd<<endl;
}
(*cord)[np][1]=(*cord)[np][1]+crp*xOverR;
(*cord)[np][3]=(*cord)[np][3]+crp*yOverR;
(*cord)[np][5]=(*cord)[np][5]+czd;	
}
of.close();
ztTozp(cord,324);
}
 */
/*
   void CSpaceCharge::getBeam(int i)
   {
   string str;
   stringstream ss;
   ss<<i;
   ss>>str;
   str+="BeamCord";
   ofstream ofile(str.c_str());
   for(int i=0;i<cord->size();i++)
   {
   for(int j=0;j<(*cord)[i].size();j++)
   ofile<<setw(8)<<(*cord)[i][j]<<" ";
   ofile<<endl;
   }
   ofile.close();
   }
   void CSpaceCharge::getvvder()
   {
   ofstream ofile("vvder");
   for(int i=0;i<vvder.size();i++)
   {
   for(int j=0;j<vvder[i].size();j++)
   ofile<<vvder[i][j]<<" ";
   ofile<<endl;
   }
   ofile.close();
   }
 */
void CSpaceCharge::getvvd(vector<vector<double> > vvd,string str,int iflag1,int iflag2)
  //compData iflag1=1,others ifalg1=0;source related ifalg=0,field related iflag2=1
{
  ofstream ofile(str.c_str());
  double dTmpr,dTmpz;
  if(iflag1==0)
  {
    for(int i=0;i<vvd.size();i++)
    {
      for(int j=0;j<vvd[i].size();j++)
      {  
        if(iflag2==0) {dTmpr=vdrs[i]-rHalf;dTmpz=vdzs[j]-zHalf;}
        else if(iflag2==1) {dTmpr=vdrm[i]-rHalf;dTmpz=vdzm[j]-zHalf;}
        ofile<<dTmpr<<" "<<dTmpz<<" "<<vvd[i][j]<<endl;
      }
    }
  }
  else if(iflag1==1)
  {
    for(int i=0;i<vvd.size();i++)
    {
      for(int j=0;j<(vvd[i].size())/2;j++)
      {  
        if(iflag2==0) {dTmpr=vdrs[i]-rHalf;dTmpz=vdzs[j]-zHalf;}
        else if(iflag2==1) {dTmpr=vdrm[i]-rHalf;dTmpz=vdzm[j]-zHalf;}
        ofile<<dTmpr<<" "<<dTmpz<<" "<<vvd[i][j]<<endl;
      }
    }
  }
}
/*
   void CSpaceCharge::getComp(string str1)
   {
//  ofstream ofile(str.c_str());
ofstream ofile(str1.c_str());
for(int i=0;i<compData.size();i++)
{
for(int j=0;j<compData[i].size();j=j+2)  //just out put the real part
//    for(int j=0;j<compData[i].size();j++)
{
//          if(compData[i][j]>1e-4)
ofile<<compData[i][j]<<" ";
//           else
//           ofile<<0<<" ";
}
ofile<<endl;
}
ofile.close();
}
 */
/*
   void CSpaceCharge::getVvv(vector<vector<vector<double> > > *vvv,string str,int iflag,double dTmpx,double dTmpy,double dTmpz,int m)
   {
   ofstream of(str.c_str());
   if(m==0)
   {
   for(int i=0;i<vvv->size();i++)
   for(int j=0;j<(*vvv)[i].size();j++)
   {
   if(iflag==0)  
   for(int k=0;k<(*vvv)[i][j].size();k++)
   {
   of<<vdxs[i]-xHalf+dTmpx/2.0<<" "<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[i][j][k]<<endl;
   }
   else if(iflag==1)
   for(int k=0;k<((*vvv)[i][j].size())/2;k++)
   {
   of<<vdxs[i]-xHalf+dTmpx/2.0<<" "<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[i][j][2*k]<<endl;
   }
   }
   }
   else if(m==1)   //输出y:z:vv[x==0][][]
   {
   for(int j=0;j<(*vvv)[nr/2-1].size();j++)
   {
   if(vdxs[nr-1]-xHalf+dTmpx/2.0<1e-5)
   {
   if(iflag==0)  
   for(int k=0;k<(*vvv)[nr/2-1][j].size();k++)
   of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[nr/2-1][j][k]<<endl;
   else if(iflag==1)
   for(int k=0;k<((*vvv)[nr/2-1][j].size())/2;k++)
   of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[nr/2-1][j][2*k]<<endl;
   }
   else 
   {
   if(iflag==0)  
   for(int k=0;k<(*vvv)[nr/2-1][j].size();k++)
   of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[nr/2-1][j][k]+(*vvv)[nr/2][j][k])/2.0<<endl;
   else if(iflag==1)
   for(int k=0;k<((*vvv)[nr/2-1][j].size())/2;k++)
   of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[nr/2-1][j][2*k]+(*vvv)[nr/2-1][j][k])/2.0<<endl;
   }
   }
   }
   else if(m==2)   //输出x:y:vv[x==0][][]
   {
   for(int j=0;j<(*vvv).size();j++)
   {
   if(vdzs[nz/2-1]-zHalf+dTmpz/2.0<1e-5||vdzs[2*(nz/2-1)]-zHalf+dTmpz/2.0<1e-5)     
   {
   if(iflag==0)  
   for(int k=0;k<(*vvv)[j].size();k++)
   of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<(*vvv)[j][k][(nz/2-1)]<<endl;
   else if(iflag==1)
   for(int k=0;k<(*vvv)[j].size();k++)
   of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<(*vvv)[j][k][2*(nz/2-1)]<<endl;
   }
   else 
   {
   if(iflag==0)  
   for(int k=0;k<(*vvv)[j].size();k++)
   of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<((*vvv)[j][k][nz/2-1]+(*vvv)[j][k][nz/2])/2.0<<endl;
   else if(iflag==1)
   for(int k=0;k<(*vvv)[j].size();k++)
   of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<((*vvv)[j][k][2*(nz/2-1)]+(*vvv)[j][k][2*(nz/2)])/2.0<<endl;
   }
   }
   }
   else if(m==3)   //输出x:z:vv[x==0][][]
   {
for(int j=0;j<(*vvv).size();j++)
{
  if(vdys[ny/2-1]-yHalf+dTmpy/2.0<1e-5)
  {
    if(iflag==0)  
      for(int k=0;k<(*vvv)[j][ny/2-1].size();k++)
        of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[j][ny/2-1][k]<<endl;
    else if(iflag==1)
      for(int k=0;k<((*vvv)[j][ny/2-1].size())/2;k++)
        of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[j][ny/2-1][2*k]<<endl;
  }
  else 
  {
    if(iflag==0)  
      for(int k=0;k<(*vvv)[j][ny/2-1].size();k++)
        of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[j][ny/2-1][k]+(*vvv)[j][ny/2][k])/2.0<<endl;
    else if(iflag==1)
      for(int k=0;k<((*vvv)[j][ny/2-1].size())/2;k++)
        of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[j][ny/2-1][2*k]+(*vvv)[j][ny/2][2*k])/2.0<<endl;
  }
}
}
of.close();
}
  */
void CSpaceCharge::getPotential(string str)
{
  ofstream of(str.c_str());
  for(int i=0;i<compData.size();i++)
  {
    for(int j=0;j<(compData[i].size())/2;j++)
      of<<vdrs[i]-rHalf<<" "<<vdzs[j]-zHalf<<" "<<compData[i][2*j]<<endl;;
  }
  of.close();
}
